SPECTRO_PEAKS

Overview

The SPECTRO_PEAKS function fits common spectroscopic peak models to experimental data, making it a versatile tool for analyzing spectral line shapes in applications such as chromatography, mass spectrometry, X-ray diffraction, and optical spectroscopy. The function uses non-linear least squares optimization via SciPy’s curve_fit to determine optimal model parameters.

Seven peak models are available, each suited for different physical phenomena:

Gaussian Peak: Models peaks where broadening arises from random processes (e.g., Doppler broadening or instrumental resolution). The peak shape follows the normal distribution.
Lorentzian Peak: Describes peaks from natural linewidth or lifetime broadening, characterized by broader tails than Gaussian profiles. Also known as the Cauchy distribution.
Gaussian-Lorentzian Dual Peak: Combines both Gaussian and Lorentzian components at the same center position, useful when both broadening mechanisms contribute to a single peak.
Pseudo-Voigt (Shared Width): A computationally efficient approximation to the true Voigt profile, expressed as a linear combination of Gaussian and Lorentzian profiles sharing a common width parameter.
Pseudo-Voigt (Independent Widths): An extended pseudo-Voigt model allowing separate Gaussian (w_G) and Lorentzian (w_L) width parameters for greater flexibility.
Voigt Profile (Convolution): The true Voigt profile computed as the convolution of Gaussian and Lorentzian distributions using the Faddeeva function. This is the most physically accurate model for spectral lines affected by both Doppler and natural broadening.
Pearson Type VII: A flexible symmetric peak model that can interpolate between Gaussian (as the shape parameter \mu \to \infty) and Lorentzian (when \mu = 1) profiles.
Inverse Polynomial Peak: Models asymmetric or complex peak shapes using a rational function with polynomial terms.

The Voigt profile is defined mathematically as:

V(x; \sigma, \gamma) = \int_{-\infty}^{\infty} G(x'; \sigma) \cdot L(x - x'; \gamma) \, dx' = \frac{\text{Re}[w(z)]}{\sigma\sqrt{2\pi}}

where w(z) is the Faddeeva function evaluated at z = (x + i\gamma) / (\sigma\sqrt{2}), \sigma is the Gaussian width, and \gamma is the Lorentzian width.

The function returns fitted parameter values along with standard errors derived from the covariance matrix, enabling uncertainty quantification for the estimated peak characteristics. For more details on the underlying optimization algorithm, see the SciPy curve_fit documentation.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=SPECTRO_PEAKS(xdata, ydata, spectro_peaks_model)

xdata (list[list], required): The xdata value
ydata (list[list], required): The ydata value
spectro_peaks_model (str, required): The spectro_peaks_model value

Returns (list[list]): 2D list [param_names, fitted_values, std_errors], or error string.

Examples

Example 1: Demo case 1

Inputs:

spectro_peaks_model	xdata	ydata
gaussian_lorentzian_dual_peak	0.01	0.18436343092326107
	2.0075	3.2739956888036987
	4.005	0.13741205361162376
	6.0024999999999995	0.01
	8	0.037490212013517314

Excel formula:

=SPECTRO_PEAKS("gaussian_lorentzian_dual_peak", {0.01;2.0075;4.005;6.0024999999999995;8}, {0.18436343092326107;3.2739956888036987;0.13741205361162376;0.01;0.037490212013517314})

Expected output:

y0	xc	A1	A2	w1	w2
0.06822	1.965	7.317	-1.657	1.706	4.96

Example 2: Demo case 2

Inputs:

spectro_peaks_model	xdata	ydata
lorentzian_peak	0.1	2.7494029905519506
	1.3250000000000002	1.0540770284615846
	2.5500000000000003	0.43012254911108183
	3.7750000000000004	0.27143948030788534
	5	0.10478140054086485

Excel formula:

=SPECTRO_PEAKS("lorentzian_peak", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.7494029905519506;1.0540770284615846;0.43012254911108183;0.27143948030788534;0.10478140054086485})

Expected output:

A	x0	w
3.022	-0.2667	1.164
0.3141	0.2236	0.06378

Example 3: Demo case 3

Inputs:

spectro_peaks_model	xdata	ydata
pearson_type_vii_peak	0.1	0.04913579064648669
	1.3250000000000002	1.7688143080814422
	2.5500000000000003	0.650921022003965
	3.7750000000000004	0.01
	5	0.013963187192617028

Excel formula:

=SPECTRO_PEAKS("pearson_type_vii_peak", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {0.04913579064648669;1.7688143080814422;0.650921022003965;0.01;0.013963187192617028})

Expected output:

xc	A	w	mu
1.72	2.463	1.162	5.417
0.07173	0.7218	0.2927	7.334

Example 4: Demo case 4

Inputs:

spectro_peaks_model	xdata	ydata
pseudo_voigt_shared_width	0.01	0.7217206195343803
	2.0075	-1.596376566988804
	4.005	0.46671030759063314
	6.0024999999999995	0.09353111787881777
	8	0.07594589700449908

Excel formula:

=SPECTRO_PEAKS("pseudo_voigt_shared_width", {0.01;2.0075;4.005;6.0024999999999995;8}, {0.7217206195343803;-1.596376566988804;0.46671030759063314;0.09353111787881777;0.07594589700449908})

Expected output:

y0	xc	A	w	mu
-0.0786	8	0.001929	0.01905	-1.935

Example 5: Demo case 5

Inputs:

spectro_peaks_model	xdata	ydata
pseudo_voigt_independent_widths	0.01	0.7217206195343803
	2.0075	-1.596376566988804
	4.005	0.46671030759063314
	6.0024999999999995	0.09353111787881777
	8	0.07594589700449908

Excel formula:

=SPECTRO_PEAKS("pseudo_voigt_independent_widths", {0.01;2.0075;4.005;6.0024999999999995;8}, {0.7217206195343803;-1.596376566988804;0.46671030759063314;0.09353111787881777;0.07594589700449908})

Expected output:

y0	xc	A	wG	wL	mu
-0.4466	5.719	31.19	0.2117	0.01549	29.75

Example 6: Demo case 6

Inputs:

spectro_peaks_model	xdata	ydata
inverse_polynomial_peak	0.5	0.27896815
	1.3	0.7237742
	2.1	2.43656316
	2.9	2.8919391
	3.7	1.01852745
	4.5	0.32536436

Excel formula:

=SPECTRO_PEAKS("inverse_polynomial_peak", {0.5;1.3;2.1;2.9;3.7;4.5}, {0.27896815;0.7237742;2.43656316;2.8919391;1.01852745;0.32536436})

Expected output:

y0	xc	w	A	A1	A2	A3
0.3058	2.4997	0.5414	0.0335	-0.5188	0.0317	-0.0004

Example 7: Demo case 7

Inputs:

spectro_peaks_model	xdata	ydata
voigt_profile_convolution	0.1	1
	1.3250000000000002	1
	2.5500000000000003	1
	3.7750000000000004	1
	5	1

Excel formula:

=SPECTRO_PEAKS("voigt_profile_convolution", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {1;1;1;1;1})

Expected output:

y0	xc	A	wG	wL
1	0.1	0	0.8167	0.4083

Python Code

import numpy as np
from scipy.optimize import curve_fit as scipy_curve_fit
from scipy.special import wofz
import math

def spectro_peaks(xdata, ydata, spectro_peaks_model):
    """
    Fits spectro_peaks models to data using scipy.optimize.curve_fit. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html for details.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        xdata (list[list]): The xdata value
        ydata (list[list]): The ydata value
        spectro_peaks_model (str): The spectro_peaks_model value Valid options: Gaussian Lorentzian Dual Peak, Lorentzian Peak, Pearson Type Vii Peak, Pseudo Voigt Shared Width, Pseudo Voigt Independent Widths, Inverse Polynomial Peak, Voigt Profile Convolution.

    Returns:
        list[list]: 2D list [param_names, fitted_values, std_errors], or error string.
    """
    def _validate_data(xdata, ydata):
        """Validate and convert both xdata and ydata to numpy arrays."""
        for name, arg in [("xdata", xdata), ("ydata", ydata)]:
            if not isinstance(arg, list) or len(arg) < 2:
                raise ValueError(f"{name}: must be a 2D list with at least two rows")
            vals = []
            for i, row in enumerate(arg):
                if not isinstance(row, list) or len(row) == 0:
                    raise ValueError(f"{name} row {i}: must be a non-empty list")
                try:
                    vals.append(float(row[0]))
                except Exception:
                    raise ValueError(f"{name} row {i}: non-numeric value")
            if name == "xdata":
                x_arr = np.asarray(vals, dtype=np.float64)
            else:
                y_arr = np.asarray(vals, dtype=np.float64)

        if x_arr.shape[0] != y_arr.shape[0]:
            raise ValueError("xdata and ydata must have the same number of rows")
        return x_arr, y_arr

    # Model definitions dictionary
    models = {
        'gaussian_lorentzian_dual_peak': {
            'params': ['y0', 'xc', 'A1', 'A2', 'w1', 'w2'],
            'model': lambda x, y0, xc, A1, A2, w1, w2: y0 + (A1 / (w1 * np.sqrt(np.pi / 2.0))) * np.exp(-2.0 * np.square((x - xc) / w1)) + (2.0 * A2 / np.pi) * (w2 / (4.0 * np.square(x - xc) + np.square(w2))),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), float(np.ptp(ya) if np.ptp(ya) else 0.5), float(np.ptp(ya) if np.ptp(ya) else 0.5), float(max((np.max(xa) - np.min(xa)) / 6.0, 1e-3)), float(max((np.max(xa) - np.min(xa)) / 6.0, 1e-3))),
            'bounds': ([-np.inf, -np.inf, -np.inf, -np.inf, 0.0, 0.0], np.inf),
        },
        'lorentzian_peak': {
            'params': ['A', 'x0', 'w'],
            'model': lambda x, A, x0, w: A * (w * w) / (np.power(x - x0, 2) + w * w),
            'guess': lambda xa, ya: (float(max(ya)), float(np.median(xa)), 1.0),
            'bounds': ([-np.inf, -np.inf, 0.0], np.inf),
        },
        'pearson_type_vii_peak': {
            'params': ['xc', 'A', 'w', 'mu'],
            'model': lambda x, xc, A, w, mu: A * np.power(1.0 + 4.0 * (np.power(2.0, 1.0 / mu) - 1.0) * np.square(x - xc) / np.square(w), -mu),
            'guess': lambda xa, ya: (float(np.median(xa)), float(np.ptp(ya) if np.ptp(ya) else 1.0), float(max((np.max(xa) - np.min(xa)) / 4.0, 1e-3)), 1.5),
            'bounds': ([-np.inf, 0.0, 0.0, 0.0], np.inf),
        },
        'pseudo_voigt_shared_width': {
            'params': ['y0', 'xc', 'A', 'w', 'mu'],
            'model': lambda x, y0, xc, A, w, mu: y0 + A * (mu * (2.0 / np.pi) * (w / (4.0 * np.square(x - xc) + np.square(w))) + (1.0 - mu) * (np.sqrt(4.0 * np.log(2.0)) / (np.sqrt(np.pi) * w)) * np.exp(-4.0 * np.log(2.0) * np.square(x - xc) / np.square(w))),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), float(np.ptp(ya) if np.ptp(ya) else 1.0), float(max((np.max(xa) - np.min(xa)) / 4.0, 1e-3)), 0.5),
            'bounds': ([-np.inf, -np.inf, -np.inf, 0.0, -np.inf], np.inf),
        },
        'pseudo_voigt_independent_widths': {
            'params': ['y0', 'xc', 'A', 'wG', 'wL', 'mu'],
            'model': lambda x, y0, xc, A, wG, wL, mu: y0 + A * (mu * (2.0 / np.pi) * (wL / (4.0 * np.square(x - xc) + np.square(wL))) + (1.0 - mu) * (np.sqrt(4.0 * np.log(2.0)) / (np.sqrt(np.pi) * wG)) * np.exp(-4.0 * np.log(2.0) * np.square(x - xc) / np.square(wG))),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), float(np.ptp(ya) if np.ptp(ya) else 1.0), float(max((np.max(xa) - np.min(xa)) / 4.0, 1e-3)), float(max((np.max(xa) - np.min(xa)) / 6.0, 1e-3)), 0.5),
            'bounds': ([-np.inf, -np.inf, -np.inf, 0.0, 0.0, -np.inf], np.inf),
        },
        'inverse_polynomial_peak': {
            'params': ['y0', 'xc', 'w', 'A', 'A1', 'A2', 'A3'],
            'model': lambda x, y0, xc, w, A, A1, A2, A3: y0 + A / (1.0 + A1 * np.power(2.0 * (x - xc) / w, 2) + A2 * np.power(2.0 * (x - xc) / w, 4) + A3 * np.power(2.0 * (x - xc) / w, 6)),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), float(max((np.max(xa) - np.min(xa)) / 6.0, 1e-3)), float(np.ptp(ya) if np.ptp(ya) else 1.0), 0.0, 0.0, 0.0),
            'bounds': ([-np.inf, -np.inf, 0.0, -np.inf, -np.inf, -np.inf, -np.inf], np.inf),
        },
        'voigt_profile_convolution': {
            'params': ['y0', 'xc', 'A', 'wG', 'wL'],
            'model': lambda x, y0, xc, A, wG, wL: (np.asarray(y0, dtype=np.float64) + np.asarray(A, dtype=np.float64) * (np.real(wofz(((np.asarray(x, dtype=np.float64) - xc) + 1j * (wL / 2.0)) / ((wG / (2.0 * np.sqrt(2.0 * np.log(2.0)))) * np.sqrt(2.0)))) / ((wG / (2.0 * np.sqrt(2.0 * np.log(2.0)))) * np.sqrt(2.0 * np.pi)))),
            'guess': lambda xa, ya: (float(np.min(ya)), float(xa[np.argmax(ya)]), max(0.0, float(np.trapz(ya - np.min(ya), xa))), float((np.max(xa) - np.min(xa)) / 6.0) or 1.0, float((np.max(xa) - np.min(xa)) / 12.0) or 0.5),
            'bounds': ([-np.inf, -np.inf, -np.inf, 0.0, 0.0], np.inf),
        }
    }

    # Validate model parameter
    if spectro_peaks_model not in models:
        return f"Invalid model: {str(spectro_peaks_model)}. Valid models are: {', '.join(models.keys())}"

    model_info = models[spectro_peaks_model]

    # Validate and convert input data
    try:
        x_arr, y_arr = _validate_data(xdata, ydata)
    except ValueError as e:
        return f"Invalid input: {e}"

    # Perform curve fitting
    try:
        p0 = model_info['guess'](x_arr, y_arr)
        bounds = model_info.get('bounds', (-np.inf, np.inf))
        if bounds == (-np.inf, np.inf):
            popt, pcov = scipy_curve_fit(model_info['model'], x_arr, y_arr, p0=p0, maxfev=10000)
        else:
            popt, pcov = scipy_curve_fit(model_info['model'], x_arr, y_arr, p0=p0, bounds=bounds, maxfev=10000)

        fitted_vals = [float(v) for v in popt]
        for v in fitted_vals:
            if math.isnan(v) or math.isinf(v):
                return "Fitting produced invalid numeric values (NaN or inf)."
    except ValueError as e:
        return f"Initial guess error: {e}"
    except Exception as e:
        return f"curve_fit error: {e}"

    # Calculate standard errors
    std_errors = None
    try:
        if pcov is not None and np.isfinite(pcov).all():
            std_errors = [float(v) for v in np.sqrt(np.diag(pcov))]
    except Exception:
        pass

    return [model_info['params'], fitted_vals, std_errors] if std_errors else [model_info['params'], fitted_vals]

Overview

Excel Usage

Examples

Python Code

Online Calculator